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PREFACE 


This report documents the research conducted to update computer models for 
dynamic fluid flow simulation of the E-l Test Stand subsystems at the NASA John C. 

S tennis Space Center. This work also involved significant upgrades to the capabilities of 
the EASY/ROCETS library though the inclusion of the NIST-12 thermodynamic 
property database and development of new control system modules. This work was 
funded under contract NAS 13-564, delivery order 143. 

The authors thank Mr. Randy Holland, Mr. Larry deQuay, and Mr. Bud Nail of 
NASA and Mr. Steve Poulton of Lockheed-Martin for their support and encouragement 
during the course of this research. 
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ABSTRACT 


Computer models for dynamic fluid flow simulation are needed to predict the 
pressures, temperatures, flow rates, etc., in present and future testing operations at the 
NASA John C. Stennis Space Center (SSC). Such simulations are used in facility design, 
test scenario development, facility modification, and facility control. The ROCket 
Engine Transient Simulation (ROCETS) package, which was initially developed by Pratt 
and Whitney for NASA Marshall Space Flight Center, and EASY5x, which is a 
commercial package developed by the Boeing Co., are the two major components which 
comprise the EASY/ROCETS dynamic fluid flow simulation package which has been 
developed by Mississippi State University (MSU) personnel for use by NASA/SSC. 
Additional code has been written to handle tasks specific to ground-test facility modeling 
such as gas bottles and pressurized liquid runtanks. 

In the present research, incorporating the complete NIST- 12 property database 
into the EASY/ROCETS library has significantly enhanced the thermodynamic property 
tables. Previously models for the low-pressure and high-pressure run systems for the E-l 
Test Stand were provided. These models have been updated based on the information 
provided by Mr. Larry deQuay. In addition, new modules have been developed to model 
the Allen-Bradley PID programmable controllers and for a valve with a time delay input. 

Items for continued development are outlined in the report. 
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SECTION 1 


INTRODUCTION 

Computer models for dynamic fluid flow simulation are needed to predict the 
pressures, temperatures, flow rates, etc., in present and future testing operations at the 
NASA John C. Stennis Space Center (SSC). Such simulations are used in facility design, 
test scenario development, facility modification, and facility control. The ROCket 
Engine Transient Simulation (ROCETS) package, which was initially developed by Pratt 
and Whitney for NASA Marshall Space Flight Center, and EASY5x, which is a 
commercial package developed by the Boeing Co., are the two major components which 
comprise the EAS Y/ROCETS dynamic fluid flow simulation package which has been 
developed by Mississippi State University (MSU) personnel for use by NASA/SSC. 
Additional code has been written (Taylor and Follett, 1994, and Follett and Taylor, 1996) 
to handle tasks specific to ground-test facility modeling such as gas bottles and 
pressurized liquid runtanks. 

In the present research, incorporating the complete NIST- 12 property database 
into the EAS Y/ROCETS library has significantly enhanced the thermodynamic property 
tables. These revisions are discussed in Section 2. 

Previously models for the low-pressure and high-pressure run systems for the E-l 
Test Stand (Follett and Taylor, 1996) were provided. These models have been updated 
based on the information provided by Mr. Larry deQuay (1997). These new models are 
presented in Section 3. 

In addition new modules have been developed to model the Allen-Bradley PID 
programmable controllers and for a valve with a time delay input. Section 4 presents the 
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development of theses new modules. Conclusions and recommendations for future work 
are presented in Section 5. Finally the notes for model development are presented in the 
Appendix. 
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SECTION 2 


THERMODYNAMIC PROPERTY TABLE 
ENHANCEMENT AND DEVELOPMENT 

The thermodynamic property tables that were programmed with the original * 
ROCETS package were meant for rocket engine simulation and, hence, had abridged 
versions of the property tables for N2 and He. Furthermore, no properties were supplied 
for Air. As need has arisen, property tables have been added in an ad hoc manner for Air 
and N2. A primary motivation for this work was the failure of the existing package with 
high pressure (14,000 psi) N2 needed to model the high-pressure 02 system on the E-l 
Test Stand (formerly know as the CTF). 

After review of the alternatives, it was decided that the best way to proceed was to 
incorporate the NIST-12 (1995) property FORTRAN codes directly into the 
EASY/ROCETS library. This was done in two stages. First, the NIST-12 routines were 
programmed into the existing EASY/ROCETS Library, “er,” for the N2 properties. This 
allows existing models to be run in the updated system without rebuilding the models and 
with only slight modification. Next the entire EASY/ROCETS library was redesigned 
from the ground up to use the NIST- 12 routines for all property computations. The 1995 
version of the NIST- 12 property routines would not handle mixed-phase states and 
performed poorly for thermodynamic states that were near the mixed-phase region — near 
the “dome.” In the course of this work, the NIST- 12 codes were modified to compute the 
thermodynamic properties of mixed-phase states. Furthermore, an error was discovered 
in the 1995 version of the FORTRAN programs. Correction of this error allows the 
revised code to compute properties near the dome with better accuracy and more reliably. 
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“er” Library Modifications 


A major problem with the original ROCETS package was the abridged nature of 
the N2 property tables. When the high-pressure 02 system on the E-l Test Stand 
(formerly know as CTF) was first modeled using EASY/ROCETS (Follett and Taylor, 
1996), the model results were inaccurate because of the limits of the N2 property tables. 
As a first step, the existing EASY/ROCETS library (“er” in the EASY/ROCETS library 
list) was modified to use the NIST- 12 property routines (corrected-expanded version 
discussed later) for the N2 property calculations. The other property tables, H2 and 02, 
were not modified. This will allow existing models to be used with high-pressure N2 
requirements; however, it is not the most efficient way to use the NIST- 12 codes and 
does not make use of all the additional capabilities that the routines make available. A 
completely redesigned library is discussed in a following section. 

Table 1 shows the modified version of the “er” library subroutine N2PROP. 
Modifications are highlighted in bold italic type. Review of the table reveals that all 
options except “P = F(RHO, H)” are computed from the NIST- 12 routines. “P = F(RHO, 
H) was omitted since it is not currently used in EASY/ROCETS. 

Figure 1 shows a test model that was used to demonstrate that the modified library 
can handle high-pressure N2. The model is very simple with a bottle blowing down 
through a valve against a backpressure in the exit module. Figure 2 shows the bottle 
pressure history during the blow down. The figure shows that the pressure history is 
smooth and continuous not stair-step as with the unmodified library (Follett and Taylor, 
1996). 
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To demonstrate the mixed-phase capabilities of the modified library, the test 
model in Figure 3 was built and executed. The inlet feeds saturated N2 vapor at 50 psia 
through the first valve into a volume that initially contains a saturated mixture of N2 
vapor and liquid at 50 psia. This volume discharges through a valve against a 
backpressure of 25 psia in the exit. Figure 4 shows the pressure and temperature history 
of the volume, “V 1 .” As the fluid in the volume discharges through the valve, the 
pressure in the volume reduces. As this pressure reduces, the inflow of saturated vapor 
from the inlet increases. Looking at the temperature plot, it is seen that the temperature 
in the volume first decreases following the saturation temperature line as the pressure 
decreases. At about 50 sec, all of the original mixed-phase fluid has been discharged or 
boiled off. After this time, the temperature increases as the fluid in the volume is 
replaced by the warmer vapor from the inlet. Figure 5 shows the computed flow rates. 

This model assumes that the mixed-phase is continuously mixed (as a mist, for 
example) not separated and that the flow of mixed phase has no effect on the flow 
properties of the valve. 

Figure 6 shows the old model (Follett and Taylor, 1996) for the high-pressure 02 
system on the E-l Test Stand. The need for high-pressure N2 in the gas pressurization 
subsystem was one of the driving forces behind the present work to upgrade the property 
computations in EASY/ROCETS. Figure 7 shows the runtank and N2 bottle pressure 
histories. The runtank pressure controller adjusts the gas control valve to maintain the 
runtank pressure at a constant level in this case. As the bottle pressure decreases, the gas 
control valve must open up more to maintain the constant ullage pressure. Figure 8 
shows the necessary control valve Cv values. Figure 9 shows a comparison of the 
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computed value of the liquid flow rate and the desired flow rate as computed by Larry 
deQuay (1996). 

Since we were modifying the “er” library, we took the opportunity to update the 
input for the runtank module. Figure 10 shows the input window. The initial conditions 
are now set by inputting the ullage vapor enthalpy and pressure, HVG and PVG, the 
liquid enthalpy, HLG, and the vapor volume of the ullage, VVG, in the “Inputs” list on 
the left of the window. This change makes the runtank input consistent with the other 
modules in the “er” library. All initial inputs come as pressure and enthalpy. 
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' Table 1. Modified N2 Property Subroutine for er Library Using NIST- 12 Routines. 


C DATA SET MSFCN2 PROP AT LEVEL 001 AS OF 08/09/90 

SUBROUTINE N2 PROP (IUPDAT , MODLOC , OPTION , IDSC , XVAR , 

* YVAR , OVAR , PNTX , PNTY , IDSCR ) 

c*********************************************************************** 

C %BEGIN CLASS N2 PROP * 

C * 

C SUBPROGRAM N2 PROP UNCLASSIFIED SID: E950 * 

C * 

C UNITED TECHNOLOGIES CORPORATION * 

C PRATT & WHITNEY * 

C WEST PALM BEACH, FLORIDA * 

C %END CLASS N2 PROP * 

Q* ********************************************************************** 

C %BEGIN PURPOSE N2 PROP * 

C * 

C THIS SUBROUTINE ACCESSES NITROGEN PROPERTIES * 

C VIA NIST_12 NITROGEN PROPERTY Programs. * 

C * 

C %END PURPOSE N2 PROP * 



%BEGIN COMMENTS N2 PROP 


THE OPTION INPUT IS INTENDED TO REPRESENT THE FUNCTIONAL 
RELATIONSHIP DESIRED. THE OPTION CODE, INDEPENDENT 
THERMODYNAMIC STATES, DEPENDENT STATE, AND MAP SUBROUTINE 
CALLED ARE SHOWN IN THE FOLLOWING TABLE. NOTE THAT THE 


7 



Table 1. Modified N2 Property Subroutine for er Library Using NIST-12 Routines. 


FORMAT OF THE OPTION ALWAYS HAS THE 
AS THE SECOND INDEPENDENT PARAMETER. 


MAP FAMILY PARAMETER 


1 

OPTION 

1 

XVAR 

1 

YVAR 

| OVAR 

| MAP 

1* 

* 

1 

P=F ( RHO , U ) 

1 

DENSITY 

1 

INTER ENG 

| PRESSURE 

ITERAT 

1* 

1 

P=F (RHO, H) 

1 

DENSITY 

1 

ENTHALPY 

| PRESSURE 

NRHP05 

1* 

1 

RHO=F (P, H) 

1 

PRESSURE 

1 

ENTHALPY 

j DENSITY 

NRHP05 

1* 

1 

T=F (P, H) 

1 

PRESSURE 

1 

ENTHALPY 

| TEMPERATURE 

NPHT05 

1* 

1 

S=F (H, P) 

1 

ENTHALPY 

1 

PRESSURE 

| ENTROPY 

NHPS05 

1* 

1 

H=F(S,P) 

1 

ENTROPY 

1 

PRESSURE 

j ENTHALPY 

NHPS05 

1* 

1 

CP=F(H,P) 

1 

ENTHALPY 

1 

PRESSURE 

j SPEC HEAT 

FLUIDS 

1* 

1 

CV=F (H, P) 


ENTHALPY 

1 

PRESSURE 

j SPEC HEAT 

FLUIDS 

1* 

1 

GAMA=F ( H , P ) 

1 

ENTHALPY 

1 

PRESSURE 

j RATIO SP 

FLUIDS 

1* 

1 

MU=F (H, P) 

1 

ENTHALPY 

1 

PRESSURE 

j VISCOSITY 

FLUIDS 

1* 

1 

K=F(H,P) 

1 

ENTHALPY 

1 

PRESSURE 

j THERMAL COND 

FLUIDS 

1* 

.__* 


%END COMMENTS N2PROP * 

**************************************** ******************************^ 

%BEGIN INTERFACE N2 PROP * 


| CALL LIST| SYSTEM 

| SYSTEM TAG 

1 

ARRAY 

| I/O 

VAR | 

| NAME 

| NAME 

1 

1 

STATUS 

j STATUS 

TYPE j 

| IUPDAT 

| IUPDAT 

| GLOBAL 

1 


| IN 0 

1*4 | 

j MODLOC 

j MODL 

j NAME 0 

1 


| IN 0 

C*4 | 

j OPTION 

j OPTN 

j VARIABLE 

1 


| IN 0 

C*12 j 

j IDSC 

j DSC 

j DSC 

1 


j IN 0 

1*4 | 

| XVAR 

j XVAR 

j VARIABLE 

1 


j IN 0 

R*4 j 

j YVAR 

j YVAR 

j VARIABLE 

1 


j IN 0 

R*4 j 

j OVAR 

j OVAR 

j VARIABLE 

1 


j OUT 0 j 

R*4 j 

j PNTX 

| PNTX 

j VARIABLE 

1 


j OUT 0 | 

R*4 j 

j PNTY 

j PNTY 

j VARIABLE 

1 


| OUT 0 j 

R*4 | 

j IDSCR 

j DSCR 

j DSCR 

1 


j OUT 0 j 

1*4 j 


C %END INTERFACE H2PROP * 

Q* *********************** *********************************************** 

C %BEGIN SUBROUTINES REQUIRED N2PROP * 

C * 

C SUBROUTINES REQUIRED : NPHT05 , NRHP05 , NHPS05 , ERCKOO, NIST_12 * 

C * 

C %END SUBROUTINES REQUIRED N2PROP * 


C %BEGIN COMMONS REQUIRED N2 PROP * 

C * 

C COMMONS REQUIRED: GUNITS, * 

C * 

C %END COMMONS REQUIRED N2 PROP * 

Q*********************************************************************** 

implicit real*8 (a-h, o-z) 

CHARACTER * 1 2 OPTION 
CHARACTER * 4 MODLOC 

COMMON / GUNITS / IUNIT , GC , GR , RJ , RU 

* CLEN , CMASS , CFORCE , CTEMP , CENERGY, 

* FLOCON 

Q* ********************************************************************* 

IF (IUNIT. EQ.O) THEN 
CFOOT=12. 
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Table 1. Modified N2 Property Subroutine for er Library Using NIST- 12 Routines. 



c 

IF ( INDEX (OPTION, 'P=F(RHO,U)' ) .NE. 0 ) THEN 

D = XVAR *CFOOT* *3 


U = YVAR 
IFL = 2 
IOP = 5 

CALL FPROP (IFL, IOP, IUNIT, P, T, QL, D, U, H, S, CP, CV, V, TB, SO, DI) 

OVAR = P 

C 

£********************************************************************** 

C 

ELSEIF ( INDEX (OPTION, 'P=F(RHO,H)' ) .NE. 0 ) THEN 

IOPT = 1 
RHO = XVAR 
H = YVAR 

CALL NRHP05 (IUPDAT , MODLOC , IOPT , IDSC , RHO 
* H , P PNTX , PNTY , IDSCR ) 

OVAR = P 

********************************************************************** 

ELSEIF ( INDEX ( OPTION , 'RHO=F(P,H)' ) .NE. 0 ) THEN 

IFL = 2 
IOP = 6 
H = YVAR 
P = XVAR 

CALL FPROP (IFL, IOP, IUNIT, P, T, QL, D, U, H, S, CP, CV, V, TH, SO, DI) 

OVAR = D/(CFOOT**3) 

********************************************************************** 

ELSEIF ( INDEX ( OPTION , 'S=F(H,P)' ) .NE. 0 ) THEN 

XFL = 2 
IOP = 6 
H = XVAR 

P = YVAR 

CALL FPROP (IFL, IOP, IUNIT, P, T, QL, D, U, H, S, CP, CV, V, TH, SO, DI) 

OVAR = S 

C 

c********************************************************************** 

c 

ELSEIF ( INDEX ( OPTION , 'H=F(S,P)' ) .NE. 0) THEN 

IOP = 4 
IFL = 2 
S = XVAR 

P = YVAR 

CALL FPROP (IFL, IOP, IUNIT, P, T, D, U, H, S, CP, CV, V, TH, SO, DI) 

OVAR = H 

C 

c********************************************************************** 

c 

ELSEIF ( INDEX ( OPTION , * T=F ( P , H ) ' ) -NE. 0 ) THEN 
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Table 1. Modified N2 Property Subroutine for er Library Using NIST- 12 Routines. 


ifl = 2 
top = 6 

P = XVAR 

H = YVAR 

CALL FPROP (IFL, TOP, IUNIT, P, T, QL, D, U, H, S, CP, CV, V, TH, SO, DI) 

OVAR = T 

C 

c** ************************************************************ ******** 

c 

ELSEIF ( INDEX ( OPTION , '=F(H,P)' ) .NE. 0) THEN 

IFL = 2 
IOP = 6 
P = YVAR 
H - XVAR 

CALL FPROP (IFL, IOP , IUNIT, P, T, QL, D, U, H, S, CP, CV, V, TH, SO, DI) 

C 

£********************************************************************** 

c 

IF (INDEX (OPTION, 'CP=').NE.O) THEN 
OVAR = CP 
C 

C********************************************** ************************ 

c 

ELSE IF (INDEX {OPTION, , CV= , ).NE.O) THEN 
OVAR = CV 
C 

Q********************************************************************** 

c 

ELSE IF (INDEX (OPTION, / GAHA= ' ) .NE . 0 ) THEN 
GAMA = CP/CV 
OVAR = GAMA 
C 

Q********************************************************************** 

C 

ELSE IF (INDEX (OPTION, 'MU=')-NE.O) THEN 
RMU = V/CVIS 
OVAR = RMU 
C 

£******************* *************************************************** 

c 

ELSE IF ( INDEX ( OPTION, ' K = 9 ) . NE . 0 ) THEN 
RK = TH/CTHK 
OVAR = RK 
C 

Q* ********************************************************************* 

C 

ELSE 

CALL ERCK00 (IUPDAT , 'N2PROP ', MODLOC//' 1ST', 10000 , 

* ' INVALID OPTION IN N2 PROP ') 

RETURN 
END IF 

ELSE 

CALL ERCK00 (IUPDAT , 'N2PROP ', MODLOC//' 2nd', 10000 

* ' INVALID OPTION IN N2 PROP ') 

GOTO 99 

ENDIF 

CALL ERCK00 (IUPDAT , 'N2PROP ', MODLOC//' 3rd', 0000 

* ' VALID OPTION IN N2 PROP ') 

C 

1000 CALL ERCK00 (IUPDAT , 'N2PROP ', MODLOC//' 4th', 0000 

* ' PRESSURE NOT FOUND IN N2 PROP — - ' ) 
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Table 1. Modified N2 Property Subroutine for er Library Using NIST- 12 Routines. 

GOTO 99 ~ 

2000 CALL ERCK00 (IUPDAT , 'N2PROP MODLOC//' 5th', 0000 

* ' ITERATIONS DONOT CONVERGE IN N2 PROP ') 

GOTO 99 

99 CONTINUE 
RETURN 
END 


C 


DATA SET MSFCERCK00 AT LEVEL 001 AS OF 08/09/90 




Fig. 1. High-Pressure N2 Property Test Model 
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Fig. 3. Mixed-Phase Test Model 
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Fig. 5. Mixed-Phase Test Model Flow Rates. 
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Fig. 10. Input Window for Revised er — Library Runtank Module. 





NIST- 12 Code Modifications 


The 1995 version of the NIST-12 property routines, which was supplied to us by 
Mr. Steve Poulton of Lockheed SSC, would not handle mixed liquid-vapor states. 
Furthermore, it would not compute saturation properties accurately in some cases. The 
codes have been modified to allow mixed-phase states. More importantly, a 
typographical error in the NIST- 12 FORTRAN program was discovered, and once 
corrected, the code now computes properties on or near the dome accurately and robustly. 

Table 2 shows the FORTRAN code for the modified subroutine FPROP. This 
subroutine drives all of the property calculations and is the one, which is directly called 
by EASY/ROCETS. The modifications are highlighted in bold italic type. The mixed 
phase calculations are based around the concept of vapor quality 

, mass vapor 

qual = 

total mass 

Therefore, qual = 1 for a saturated vapor, and qual =0 for a saturated liquid. For mixed 
states 0 < qual < 1.0. Here the convention is established that for all nonsaturated states 
such as gases, compressed liquids, or supercritical states qual will be set to -1 (qual = - 
1.0). In the mixed states all intensive properties (specific volume, 1/p, internal energy, u, 
enthalpy, h, and entropy, s) are computed using the quality 

— = qual h (1 — qual) 

P P v Pl 

u = qual • u v + (1 - qual) • u L 
h = qual -h v +( 1 - qual) ■ h L 
s = qual • s v + (1 - qual) • s L 
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This is not necessarily true for the mixed phase values of the other properties. Specific 
heats, c p and c v , and viscosity may not have meaning for mixed phases. However, they 
are naively computed using the formulas above. The mixed-phase values for di, cp, cv, v, 
th, and so should not be used without further study. 

When calls are made for saturated properties at a given pressure or temperature 
the user must supply the quality. When calls are made using intensive properties, such as 
pressure and internal energy or density and internal energy, the code tests for saturated 
conditions and computes quality. 

Very small modifications to the FORTRAN functions FINDP and FINDT were 
also necessary and are shown in Tables 3 and 4. 

Table 5 shows the FORTRAN code for the function VPN. This function 
computes the vapor pressure at a given temperature. The highlighted line, shown here in 
the corrected form, had a typographical error. The term VP(2) was repeated twice as 
VP(2)*VP(2). This caused significant inaccuracies for some materials, e.g. hydrogen. 

The uncorrected code would not compute saturation properties correctly and was not 
robust computing properties near the dome since the improper saturation properties were 
being used as starting guesses in many of the iterative functions, FINDD for example. 

Table 6 shows before and after results for saturated hydrogen vapor at 50 psia 
along with values from the “Hydrogen Technological Survey Thermophysical Properties” 
tables (McCarty, 1975). From the table, we see that the uncorrected routines failed badly. 
First the saturation temperature is computed incorrectly. This error causes mass 
confusion with the other routines and the properties come out to be those for compressed 
liquid at 50.0 psia and 35.6 R. The computed properties after the correction agree very 
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closely with those given in the Tables. The energy values, u and h, are in slight error 
because the routines compute these saturated values at very small superheat to maintain 
the robustness of the iterative procedures imbedded in the code. The old “miprops” 
program does the same thing and gets very close to the same results as the after column 
in the table. 

Figure 1 1 shows the input stream of the NIST- 12 program for a mixed-phase state 
of N2 at 50 psia and a quality of 0.5. Figure 12 shows the results for the same state when 
the inputs are density and internal energy. Notice that in this case a starting temperature 
guess is now required. If this option has trouble converging, try a different temperature 
guess. 

The corrected code has proven to be very robust in that it usually converges. The 
one exception is near the critical point. Figure 13 shows a temperature-density-pressure 
map for parahydrogen (McCarty, 1975). The heart of the NIST database is a curve fit of 
this map in the form p = f(t, p). This formula is used in an iterative mode in several 
places in the code. Inspections of the figure show that the isobars have an inflection 
point at the critical point. This causes the iterative algorithms to be unstable very near 
the critical point. Two fixes have been thought of. One is to box the critical point with 
four nearby points and to use interpolation when inside the box. The other is to use a bi- 
section iterative scheme instead of the Newton’s method currently used in NIST- 12. 
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Table 2. FORTRAN Code for Modified NIST 12 FPROP Subroutine 


£**************************★** Subiroxi t ine ****************************** 
subroutine fprop (ifl, iop, iu,p, t,qual,d,u,h, s,cp,cv,v, 

& th,so,di) 

implicit real*8 (a-h, o-z) , integer*4 ( i-n) 


common 

/crit/ 

em, 

eok. 

rm, 

tc. 

dc. 

Sc 


x. 

pc. 

sig 



common 

/data/ 

g(32). 

r. 

gam, 

vp(9) , 

dtp, 

& 


pec, 

Ptp, 

tcc. 

ttp. 

tul. 

& 


til. 

pul. 

dec 



common 

/con/ 

peon. 

tcon. 

dcon. 

ucon, 

scon 

& 


vc on. 

thcon. 

socon 



common 

/satf 1/ 

il, 

ip. 

isat. 

dv. 

dl 


End of COMMON Structures 


loop = 2 

if (ifl.eq.l) then 
call PH2 ( ) 

else if (ifl.eq.2) then 
call N2{) 

else if (ifl.eq.3) then 
call 02() 

else if (ifl.eq.4) then 
call AR() 

else if (ifl.eq.5) then 
call NF3 () 

else if (ifl.eq.6) then 
call CH4 ( ) 

else if (ifl.eq.7) then 
call C2H6 0 

else if (ifl.eq.8) then 
call C2H4() 

else if (ifl.eq.9) then 
call C3H8() 

else if (ifl.eq.10) then 
call ISOB() 

else if (ifl.eq.il) then 
call NORBO 

else if (if 1 . eq. 12 ) then 
call D2() 

else if (ifl.eq.13) then 
call HE ( ) 

else if (ifl.eq.14) then 
call C02 ( ) 

else if (ifl.eq.15) then 
call COO 

else if (ifl.eq.16) then 
call NH2() 

else if (ifl.eq.17) then 
call XE ( ) 

else 

write ( * , * ) ' UNKNOWN FL 

return 

end if 


' UNKNOWN FLUID NUMBER PLEASE REENTER' 


C + + H 

if (iu . 

eq . 0 ) then 


c 

Engr . 

Units - Conversions 

* 


peon 

= 14.695949d0/0. 101325d0 



tcon 

= 1 . 8d0 



dcon 

= em/16. 0184637d0 



ucon 

= 1 . / (2 . 324445d0*em) 



scon 

= 1 . / (4 . 184001d0*em) 



vc on 

= 1 . d-6*0 . 671875652d0 
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Table 2 . FORTRAN Code for Modified NIST_12 FPROP Subroutine 


thcon = 0 . 578176d0 
socon = 3 . 28084d0 
else 

SI Units - Conversions 
peon = 1 . dO 
tcon = l.dO 
dcon = em 

ucon = 1 .dO/em*l . 0d3 
scon = 1 . dO/em*l . 0d3 
vcon = 1.0d-3 
thcon = 1 . dO 
socon = l.dO 
endif 

write (*,*) pcc*pcon, tcc*tcon 
PP = p/pcon 
tt = t/tcon 
dd = d/dcon 
uu = u/ucon 
hh - h/ucon 
ss = s/scon 
epp = cp/scon 
cw = cv/scon 
soo = so/socon 
w = v/vcon 
thh = th/ thcon 


if {iop . eq. 0) then 

call satur( pp, tt, ip ) 

if (il .It. 0) then 
plus = 0 . 025&-Q2 *pp 
dl = £indd(pp+plus,tt) 
dv findd(pp-plus,tt) 
call energy ( pp, tt , dv, uv, hv, sv ) 
call spheat ( tt, dv, eppv, cvw, soov ) 
div = dielec ( pp, tt, dv ) 
call visthe( dv, tt, wv, tbhv ) 
call energy( pp, tt, dl, ul, hi, si ) 
call spheat ( tt, dl, cppl, cwl, sool ) 
dll = dielec ( pp, tt, dl ) 
call visthe ( dl, tt, ml, thhl ) 
dd=qual * (1 . OdO/dv) + (1 . OdO-gual )*(1. OdO/dl ) 
dd=l . OdO/dd 

uu-gual *uv+(l. OdO-qual) *ul 
hh=qual *hv+ (1 . OdO-qual) *hl 
ss=qual *sv+(l . OdO-qual) *sl 
di=:qual *div+ (1 . OdO-qual) *dil 
cpp=qual*cppv+ ( 1 . OdO-qual) *cppl 
cm-qual *cvm+ (1 . OdO-qual ) *cml 
m=qual *vm+(l. OdO-qual) *ml 
thh=qual *thhv+ (1 . OdO-qual ) *thhl 
soo=qual *soov+ (1 . OdO-qual ) *sool 
else 

if (il .eq. 0) then 
plus = 0 . 25d-01*pp 
dd = f indd( pp+plus, tt ) 

qual = O.OdO 

else 

plus = 0 . 25d-02*pp 

dd = findd( pp-plus, tt ) 

qual = l.OdO 

endif 

call energy ( pp, tt, dd, uu, hh, ss ) 
call spheat ( tt, dd, epp, cw, soo ) 
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di = dielec ( pp, tt, dd ) 
call vis the ( dd, tt, w, thh ) 
end if 

else if (iop . eq. 1) then 
call limits ( pp, tt ) 
dd = findd( pp, tt ) 

call energy ( pp, tt, dd, uu, hh, ss ) 
call spheat< tt, dd, cpp, cw, soo ) 
di = dielec ( pp, tt, dd ) 
call visthe( dd, tt, w, thh ) 

qual = -1,0 

else if (iop .eq. 2) then 

il = -1 
isat = 0 

tt - findt( pp, dd ) 
call limits ( pp, tt ) 

if (isat.eq.l) then 

call energy ( pp, tt, dv, uv, hv, sv ) 
call spheat ( tt, dv, cppv, cvw, soov ) 
div = dielec ( pp, tt, dv ) 
call visthe( dv, tt, vw, thhv ) 
call energy( pp, tt, dl, ul, hi, el ) 
call spheat ( tt, dl, cppl, cwl, sool ) 
dll = dielec ( pp, tt, dl ) 
call visthe( dl, tt, wl, tbhl ) 
qual = (l.OdO/dd - l.OdO/dl) / (l.OdO/dv - l.OdO/dl) 
uu=qual *uv+ ( 1 . OdO -qual ) *ul 
hh=qual *hv+ ( 1 • OdO -qual) *hl 
ee=qual *sv+(l . OdO -qual) *el 
di-qual *div+ (1 . OdO-qual) *dil 
cpp=qual *cppv+ ( 1 , OdO-qual) *cppl 
cw=qual *cwv+(l, OdO-qual) *cwl 
w =qual *vw+ ( 1 . OdO-qual) *wl 
thhssqual *tbhv+ ( 1 . OdO-qual) *thhl 
soo -qual *soov+(l . OdO-qual) *sool 
else 

call energy ( pp, tt, dd, uu, hh, ss ) 
call spheat ( tt, dd, cpp, cw, soo ) 
di = dielec ( pp, tt, dd ) 
call visthe{ dd, tt, w, thh ) 

qual = -1. OdO 

end if 

else if (iop .eq. 3) then 

il m -1 
isat = 0 

pp = findp( tt, dd ) 
call limits ( pp, tt ) 

if (isat,eq,l) then 

call energy ( pp, tt, dv, uv, hv, sv ) 
call spheat ( tt, dv, cppv, cvw, soov ) 
div = dialed pp, tt, dv ) 
call visthe( dv, tt, vw, thhv ) 
call energy( pp, tt, dl, ul, hi, si ) 
call spheat ( tt, dl, cppl, cwl, sool ) 
dil * dielec ( pp, tt, dl ) 
call visthe( dl, tt, wl, thhl ) 
qual * (l.OdO/dd - 1. OdO/dl) / (1. OdO/dv - l.OdO/dl) 
uu=qual *uv+(l . OdO-qual) *ul 
hh=qual *hv+ (1 . OdO-qual) *hl 
ss=qual *sv+ ( 1 • OdO-qual) *sl 
di-qual *div+ (1 . OdO-qual) *dil 
cpp-qual *cppv+ (1 . OdO-qual ) *cppl 
cw=qual *cvw+ (1 . OdO-qual) *cwl 
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w-qual *vw+(l. OdO-qual) *wl 
thh=qual *t hhv+ (1 . OdO-qual ) *thhl 
soo=qual *soov+ (1 . OdO-qual ) *sool 
else 

call energy ( pp, tt, dd, uu, hh, ss ) 
call spheat( tt, dd, cpp, cw, soo ) 
di = dielec ( pp, tt, dd ) 
call visthe{ dd, tt, w, thh ) 
qual = -l.OdO 
end if 

else if (iop .eq. 4) then 

Isat = 0 

if( (pp.lt .pcc) .and. (pp.gt.ptp)) then 
ip - 1 

call satur(pp,tsat,ip) 
dv - findd(pp-1.0d-04,tsat ) 
dl = findd(pp+0.25d-02*pp,tsat) 
call energy ( pp, tsat, dv, uv, hv, sv ) 
call energy( pp, tsat, dl, ul, hi, si ) 
if ( (ss.lt. sv) .and. (ss.gt.sl) ) isat = 1 
end if 

if (isat.eq.l) then 
tt - tsat 

call spheat( tt, dv, cppv, cvw, soov ) 
div = dielec ( pp, tt, dv ) 
call visthe( dv, tt, vw, thhv ) 
call spheat( tt, dl, cppl, cwl, sool ) 
dil m dielec ( pp, tt, dl ) 
call visthe( dl, tt, wl, thhl ) 
qual = (ss - si) / (sv - si) 
uu=qual *uv+(l. OdO-qual) *ul 
hh=qual *hv+ (1 . OdO-qual) *hl 
dd=qual * (1 . OdO/dv) + (1 . OdO-qual )*(1. OdO/dl ) 
dd=l. OdO/dd 

di=qual *div+ (1 . OdO-qual ) *dil 
cpp=qual *cppv+ (1 . OdO-qual) *cppl 
cw=qual *cvw+ (1 . OdO-qual ) *cwl 
w=qual *vw+ (1 . OdO-qual) *wl 
thh=qual *thhv+ (1 . OdO-qual) *t hhl 
soo=qual*soov+(l. OdO-qual) *sool 
else 

si = ss 
isig = 0 

do while {isig .It. 4) 
dd = findd( pp, tt ) 

call energy ( pp, tt, dd, uu, hh, ss ) 
err = ss-si 

call nuiter( tt, err, np, isig, loop ) 
end do 

if (isig .eq. 4) then 

write (*,*) 'FPROP - Option ',iop,' Failed to Converge' 
endif 

dd = findd( pp, tt ) 
call limits ( pp, tt ) 
call energy ( pp, tt, dd, uu, hh, ss ) 
call spheat( tt, dd, cpp, cw, soo ) 
di = dielec ( pp, tt, dd ) 
call vis the { dd, tt, w, thh ) 
end if 

else if {iop .eq. 5) then 
ui = uu 
isig = 0 

il = - 1 


\ 
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isa t = 0 

do while {isig .It. 4) 
pp = findp( tt, dd ) 

if (isat.eq.l) then 

call energy (pp, tt, dv, uv, hv, sv) 
call energy (pp,tt,dl,ul, hi, si) 
qual = (l.OdO/dd - l.OdO/dl) / (l.OdO/dv - l.OdO/dl) 
uu-qual*uv+(l. OdO-qual) *ul 
else 

call energy ( pp, tt, dd, uu, hh, ss ) 
end if 
err = uu-ui 

call nuiter ( tt, err, np, isig, loop ) 
end do 

if (isig .eq. 4) then 

write <*, *) 'FPROP - Option ' , iop, ' Failed to Converge' 
endif 

pp = f indp ( tt, dd ) 
call limits ( pp, tt ) 

if (isat.eq.l) then 

call energy ( pp, tt, dv, uv, hv, sv ) 
call spheat ( tt, dv, cppv, cm, soov ) 
div » dielec ( pp, tt, dv ) 
call visthe( dv, tt, vw, thhv ) 
call energy( pp, tt, dl, ul, hi, si ) 
call spheat ( tt, dl, cppl, cwl, sool ) 
dil - dielec ( pp, tt, dl ) 
call visthe( dl, tt, wl, thhl ) 
qual = (l.OdO/dd - l.OdO/dl) / (l.OdO/dv - l.OdO/dl) 
uu=qual *uv+ ( 1 . OdO-qual) *ul 
hh=qual *hv+ ( 1 . OdO-qual) *hl 
ss~qual *sv+ (1 . OdO-qual) *sl 
di=qual *div+(l. OdO-qual) *dil 
cpp=qual *cppv+ (1 . OdO-qual) *cppl 
cw=qual *cvw+ ( 1 . OdO-qual) *cwl 
w=qual *wv+ (1 • OdO-qual) *wl 
thh=qual *thhv+ (1 . OdO-qual) *t hhl 
soo-qual *soov+ ( 1 • OdO-qual) *sool 
else 

call energy ( pp, tt, dd, uu, hh, ss ) 
call spheat ( tt, dd, cpp, cw, soo ) 
di « dielec ( pp, tt, dd ) 
call vis the ( dd, tt, w, thh ) 

qual « -1 . OdO 

end if 

else if (iop .eq. 6) then 

isat = 0 

if ((PP* It.pcc) .and. (pp.gt.ptp) ) then 
ip B 1 

call satur(pp,tsat,ip) 
dv a findd (pp-1 . 0d-04, tsat) 
dl = findd(pp+0.25d-02*pp, tsat) 
call energy( pp, tsat, dv, uv, hv, sv ) 
call energy( pp, tsat, dl, ul, hi, si ) 
if((hh.lt.hv).and.(hh.gt.hl)) isat = 1 
end if 

if (isat.eq.l) then 
tt a tsat 

call spheat ( tt, dv, cppv, cvw, soov ) 

div a dielec ( pp, tt, dv ) 

call visthe ( dv, tt, vw, thhv ) 

call spheat ( tt, dl, cppl, cwl, sool ) 

dil a dielec ( pp, tt, dl ) 
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call visthe( dl , tt, w 1, thhl ) 
q ual = (hh - hl)/(hv - hi) 
uu=qual *uv+ (l.OdO -qual ) *ul 
ss=qual *sv+ ( 1 • OdO-qual) * b 1 
dd=qual * (1 . OdO/dv) + (1 . OdO -qual )*(1.0d0/dl) 
dd=l . OdO/dd 

dl=qual *div+(l. OdO-qual) *dil 
cpp=qual *cppv+ (1 . OdO -qual ) *cppl 
cw=qual *cwv+ ( 1 . OdO-qual) *cwl 
w=qual *v~w+ ( 1 . OdO-qual) *wl 
thh=qual *thhv+ (1 . OdO-qual) *thhl 
soo=qual *soov+(l . OdO-qual) *sool 
else 

hi = hh 
isig = 0 

do while (isig .It. 4) 
dd = findd( pp, tt ) 

call energy ( pp, tt, dd, uu, hh, ss ) 
err = hh-hi 

call nuiter( tt, err, np, isig, loop ) 
end do 

if (isig .eq. 4) then 

write (*, *) 'FPROP - Option ',iop,' Failed to Converge' 
endif 

dd = f indd( pp, tt ) 
call limits ( pp, tt ) 
call energy ( pp, tt, dd, uu, hh, ss ) 
call spheat( tt, dd, cpp, cw, soo ) 
di = dielec ( pp, tt, dd ) 
call visthe( dd, tt, w, thh ) 
qual = - l.OdOO 
end if 

else if (iop .eq. 7) then 
si = ss 
isig = 0 

11 = -1 
isat = 0 

do while (isig .It. 4) 
pp = findp( tt, dd ) 

if (lsat.eq.1) then 

call energy (pp, tt, dv, uv, hv, sv) 
call energy (pp,tt,dl,ul, hi, si) 
qual ~ (1, OdO/dd - l.OdO/dl) / (1. OdO/dv - l.OdO/dl) 
ss=qual *sv+ (1 . OdO-qual) *sl 
else 

call energy ( pp, tt, dd, uu, hh, ss ) 
end if 
err = ss-si 

call nuiter( tt, err, np, isig, loop ) 
end do 

if (isig .eq. 4) then 

write (*, *) 'FPROP - Option ',iop,' Failed to Converge' 
endif 

isat = 0 

pp = findp( tt, dd ) 
call limits ( pp, tt ) 

if (isat.eq.l) then 

call energy( pp, tt, dv, uv, hv, sv ) 
call sphea t( tt, dv, cppv, cvw, soov ) 
div = dielec ( pp, tt, dv ) 

call visthe( dv, tt, wv, thhv ) 
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call energy ( pp, tt, dl, ul, hi, el ) 
call spheat ( tt, dl, cppl, cwl, sool ) 
dil sr dielec ( pp, tt, dl ) 
call vis the ( dl, tt, wl, thhl ) 
gual m (l.OdO/dd - l.OdO/dl) / (l.OdO/dv - l.OdO/dl) 
uu-qual *uv+ (1 . OdO-qual) *ul 
hh-qual *hv+ ( 1 . OdO-qual) *hl 
ss=qual *sv+ (l.OdO - qual ) *sl 
di=qual *div+(l. OdO-qual) *dil 
cpp=qual *cppv+(l. OdO-qual) *cppl 
cw^qual *cvw+ (1 . OdO-qual ) *cwl 
w=qual*vw+ ( 1 . OdO-qual) *wl 
thh=qual*thhv+(l. OdO-qual) *thhl 
soo-qual *soov+ (1 . OdO-qual ) *sool 
else 

call energy( pp, tt, dd, uu, hh, ss ) 
call spheat ( tt, dd, cpp, cw, soo ) 
di = dielec ( pp, tt, dd ) 
call visthe( dd, tt, w, thh ) 

qual = -1 . OdO 

end if 

else if (iop .eq. 8) then 

isat - 0 

if( (pp.lt.pcc) . and . then 

ip *= 1 

call satur(pp,tsat,ip) 
dv = £indd(pp-l. 0d-04, teat) 
dl ~ £indd(pp+0.25d-02*pp, teat) 
call energy ( pp, teat, dv, uv, hv, ev ) 
call energy ( pp, teat, dl, ul, hi, el ) 
if ( (uu.lt. uv) .and. (uu.gt.ul) ) isat » 1 
end if 

if (isat.eq.l) then 
tt ^ teat 


call 

spheat ( tt. 

dv. 

cppv. 

CVW, 

soov ) 

div ■ 

= dielec ( pp. 

tt. 

dv ) 



call 

vie the ( dv. 

tt. 

vw. 

thhv ) 


call 

spheat ( tt. 

dl. 

cppl. 

cwl. 

sool ) 

dil * 

= dielec ( pp. 

tt. 

dl ) 



call 

visthe( dl. 

tt. 

wl. 

thhl ) 



qual = (uu - ul)/(uv - ul) 

dd-qual * (1 . OdO/dv) + (1 . OdO-qual) * (1 . OdO/dl) 
dd=l . OdO/dd 

hh=qual *hv+ (1 . OdO-qual) *hl 
se=qual *ev+ (1 . OdO-qual ) *sl 
di=qual *div+ (1. OdO-qual) *dil 
cpp=qual *cppv+ (1 . OdO-qual) *cppl 
cw=qual *cvw+ (1 . OdO-qual ) *cwl 
w=qual *vw+ (1 . OdO-qual) *wl 
thh=qual *thhv+ (1 . OdO-qual) *t hhl 
soo-qual *eoov+ (1 . OdO-qual ) *sool 
else 

ui = uu 
isig = 0 

do while (isig .It. 4) 
dd = findd( pp, tt ) 

call energy ( pp, tt, dd, uu, hh, ss ) 
err = uu-ui 

call nuiter( tt, err, np, isig, loop ) 
end do 

if (isig .eq. 4) then 

write (*,*) ' FPROP - Option ',iop,' Failed to Converge 7 

endif 
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dd = 

f indd ( ; 

PP/ 

tt ) 



call 

limits ( 

PP/ 

tt 

) 


call 

energy ( 

PP / 

tt. 

dd. 

uu, hh, ss 

call 

spheat ( 

tt. 

dd, 

cpp, 

, cw, soo 

di = 

dielec ( 

PP/ 

tt. 

dd ) 

i 

call 

vis the ( 

dd, 

tt. 

w. 

thh ) 


Qual = -l.dOO 

end if 
else 

pp = 0 . dO 
tt = O.dO 
dd = 0 . dO 
uu = O.dO 
hh = 0 . dO 
ss = O.dO 
cpp = O.dO 
cw = O.dO 
w = O.dO 
thh = O.dO 
soo = O.dO 
di = 0 . dO 
endif 



\ 
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Table 4. FORTRAN Code for Modified NIST_12 FINDT Function 


c ******************************* Func t ion ****************************** 
real*8 function findt (p, d) 
implicit real*8 (a-h, o-z) , integer *4 ( i-n) 


common 

/data/ 

g(32), 

ptp, 

dec 

r , 

tcc, 

gam, 

ttp. 

vp (9) , 
tul. 

dtp, 

til. 

pcc , 
pul. 

common 

/con/ 

peon, 

thcon, 

tcon, 

socon 

dcon, 

ucon. 

scon, 

vc on 

common 

/satf 1/ 

il. 

iP/ 

isat. 

dv. 

dl 



C++++++++++++++++++++++ end of COMMON Structures +++++++++++++++++++++-* 

c* Returns Temperature (K) , from the 3 2 -term MBWR eqn of state, 
c* input is Pressure (mPa) and Density (mol /I ) . 

ieq = 1 
loop = 1 

isat - 0 

if (p .gt. pcc) then 
tt = tcc 
else 

tsat = tcc 
isig = 0 

do while (isig .It. 4) 

PP = vpn ( tsat ) 
err = p - pp 

call nuiter { tsat, err, np, isig, loop ) 
end do 

if (isig .eq. 4) then 
write (*, *) 

write (*,*) 'findt - Vapor pressure equ. Failed to converge’ 

write (*,*) ' PRESS = * ,p*pcon, 1 TSAT = ' ,tsat*tcon 

write(*,*) ' ERROR = ' , (p-pp) *pcon 

write(*, *) 
endif 

PP = P " l.d-4 
dv = findd( pp, tsat ) 
pp = p + 0.25d-02*p 
dl = findd( pp, tsat ) 

if ( (d .gt. dv) .and. (d .It. dl) ) then 
isat = 1 
write (*, *) 

write(*,*) 'findt - The state point specified corresponds' 
write(*,*) ' to a Density in the 2 -phase region' 

write (*, *) 
findt = tsat 
return 
endif 
tt = tsat 
endif 
isig = 0 

do while (isig .It. 4) 

pp =s props ( d, tt, ieq ) 
err = p - pp 

call nuiter ( tt, err, np, isig, loop ) 
end do 

if (isig .eq. 4) then 
write (*, *) 

write ( * , * ) 'findt - Equation of state Failed to converge' 
write (*,*) ' PRESS = ' ,p*pcon, 1 TEMP = 1 ,tt*tcon 
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write (*,*) ' 

ERROR = ' , (p-pp) *pcon 

write ( *, *) 


endif 


findt = tt 


return 


end 


0**************************** 

end Function ***************************** 
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Table 6. Comparison of Corrected and Uncorrected Computations for Saturated 

Hydrogen Vapor at 50 psia. 



Tables 

Before 

After 

Pressure [psia] 

50.0 

50.0 

50.0 

Temperature [R] 

45.4 

35.6 

45.4 

Density [lbm/ft J ] 

0.262 

4.47 

0.261 

Internal Energy [BTU/lbm] 

52.3 

-113.2 

52.8 

Enthalpy [BTU/lbm] 

87.6 

- 111.1 

fcS 88.2 

Entropy [BTU/lbm-R] 

6.29 

1.84 

6.30 
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NIST Standard Reference Database #12 
NIST Thermophysical Properties of Pure Fluids Database 


This Program Provides the Thermodynamic Properties of 17 Fluids 


Para Hydrogen, 
Argon , 

Ethane, 

Iso Butane, 
Helium, 

Normal Hydrogen, 


Nitrogen, 

Nitrogen Trifluoride, 
Ethylene, 

Normal Butane, 

Carbon Dioxide, 

Xenon 


Oxygen , 

Methane, 

Propane , 
Deuterium, 

Carbon Monoxide, 


* Original Development by R.D. McCarty and V. Arp * 

* Based on correlations described by NBS Technical Note 1097 * 

* * 

* Revised - May 15, 1995 * 

* * 

************************************************************************ 
Is the Input comming from a file? (y/n) 

Default - (n) : 


Select A Fluid By Entering The Corresponding Number 
************************************************************************ 

* l=Para Hydrogen 2=Nitrogen 3-Oxygen * 

* 4=Argon 5=Nitrogen Tri fluoride 6=Me thane * 

* 7 =E thane 8=Ethylene 9= Propane * 

* 10=Iso Butane ll=Normal Butane 12=Deuterium * 

* 13=Helium 14=Carbon Dioxide 15=Carbon Monoxide * 

* 16=Normal Hydrogen 17=Xenon 18=Stop * 

************************************************************************ 


Default (18) - Stop : 


For Engineering Units Enter "0 K , for Metric Units Enter M 1 M 
Default 0 - Engineering Units : 


Enter The Input Option Desired 

************************************************************************ 

* 0=Saturation Properties l=Input Pressure & Temperature * 

* 2= Input Pressure & Density 3= Input Temperature & Density * 

* 4=Input Pressure & Entropy 5=Input Density & Internal Energy * 

* 6= Input Pressure & Enthalpy 7= Input Density & Entropy * 

* 8= Input Pressure & Internal Energy * 


Default 1 - Pressure & Temperature : 

0 

For Liquid Enter M 0 H , For Vapor Enter M 1 M 
For Mixture Enter H -1 H 

Default 0 - Liquid : 

-2 

Input Temperature Enter "0", Input Pressure M l" 
Default 0 - Temperature : 

1 


To Select Another Fluid Enter Zero (0) for the Input 


Figure 1 1. Output Listing for Modified NIST_12 Property Routines with Saturated Fluid Capability- 
Pressure and Quality as Input. 


37 






Enter a Pressure in (psia) 

(Decimal Point Required) 

50 . 

Enter the Quality 

(Decimal Point Required) 

.5 




* 

Press = 

50.0 

(psia) 


* 

* 

Temp 

161.11 

(deg R) 


* 

★ 

Temp = 

-298.56 

(deg R) 


* 

★ 

Qual = 

0.500 



* 

★ 

Dens 

1.768 

(lbm/ft3) 

Z = 0.45818 

* 

* 

IntEng = 

-7.60 

(Btu/lbm) 


* 

* 

Enthal = 

-2.36 

(Btu/lbm) 


★ 

* 

Entrop = 

0.992 

(Btu/lbm-R) 


* 

★ 

Cv 

0.209 

(Btu/lbm-R) 


* 

* 

Cp 

0.404 

(Btu/lbm-R) 

k = 1.93094 

* 

* 

Sound = 

1532. 

(ft/sec) 


* 

* 

Vise = 

3.57968 

xE-5 (lbm/ft-sec) 

= 0.05328 (c-Poise) 

* 

* 

Th Con = 

0.03461 

(Btu/ft-hr-R) 


* 

* 

Diel = 

1.20256 



* 




Figure 11. Output Listing for Modified NIST_12 Property Routines with Saturated Fluid Capability- 
Pressure and Quality as Input. 
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Enter The Input Option Desired 

************************************************************************ 

* 0=Saturation Properties l=Input Pressure & Temperature * 

* 2=Input Pressure & Density 3=Input Temperature & Density * 

* 4=Input Pressure & Entropy 5=Input Density & Internal Energy * 

* 6=Input Pressure & Enthalpy 7=Input Density & Entropy * 

* 8= Input Pressure & Internal Energy ' * 

************************************************************************ 

Default 1 - Pressure & Temperature : 

5 

************************************************************************ 
To Select Another Fluid Enter Zero (0) for the Input 

Enter Density (lbm/ft3) and Internal Ergy (Btu/lbm) Temp Guess (deg R) 
(Decimal Points and Comma Required) 

1 . 768 ,- 7 . 6 , 170 . 


Press 

Temp 

Temp 

Qual 

Dens 

IntEng 

Enthal 

Entrop 

Cv 

Cp 

Sound 
Vise 
Th Con 
Diel 


50.0 

161.11 

-298.56 

0.500 

1.768 

-7.60 

-2.37 

0.992 

0.209 

0.404 

1532. 

3.57973 

0.03461 

1.20255 


(psia) 

(deg R) 

(deg R) 

(lbm/f t3 ) 

(Btu/lbm) 

(Btu/lbm) 

(Btu/lbm-R) 

(Btu/lbm-R) 

(Btu/lbm-R) 

(ft/sec) 

xE-5 (lbm/ft-sec) 
(Btu/f t-hr-R) 


= 0.45809 


= 1.93089 


= 0.05328 (c-Poise) 



Figure 12. Output Listing for Modified NIST_12 Property Routines with Saturated Fluid Capability- 
Density and Internal Energy as Input. 
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Fig. 13. Parahydrogen Temperature-Density-Pressure Map 






“nr” Library— Revised For NIST- 12 


Using the NIST- 12 routines to make all property calculations allows opportunities 
for initial condition input that were not available with the ROCETS table lookup routines. 
Also, in the original ROCETS routines each property was a separate table lookup call. 
The NIST- 12 routines are organized to compute all of the properties with each call to 
FPROP. Therefore, it would be inefficient to paste the NIST- 12 codes into the old 
ROCETS paradigm as was done in the revised “er” library discussed above. This 
prompted us to completely revise the modules to be used with the NIST- 12 routines. 

This new library is abbreviated “nr” in the EASY5x library list. 

The organization of the new library can be fully understood by reviewing four 
modules. The first is the bottle module. Figure 14 shows the input window for the bottle 
module. Looking at the “Inputs” column, it is seen that the initial state is input by 
specifying the temperature, pressure, and quality (TVG, PVG, and VQG) instead of the 
pressure and enthalpy in the “er” library. Pressure and temperature were selected as the 
more natural properties to define the initial state. When the quality is set equal to -1, the 
default, both pressure and temperature are used to compute the state properties. When 
quality is set between 0 and 1, saturation properties are computed at the specified 
pressure, and the temperature input is ignored. The fluid (H2, N2, 02) is chosen with the 
map variable. Care should be taken with the map variable. In line with the option 
numbers in NIST-12 (1=H2, 2=N2, 3=02) the map numbers have changed form the “er” 
library. Table 7 shows the EASY5 macro FORTRAN code for the nr-library bottle. 

Figure 15 shows the input window for the inlet module. Again the Input column 
contains pressure, temperature, and quality to set the initial state properties. Table 8 
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contains the FORTRAN macro code. Figure 16 shows the input window for pipe module 
PIPE99, PA. It should be noted that the Input column no longer has a place for the map 
variable since the viscosity is now fed in through the inlet port from the neighboring 
volume-type element. Table 9 contains the FORTRAN code. 

Figure 18 shows a test model for some pipe and valve elements along with the 
flow rate obtained with both the old “er” library and the new “nr” library. The results are 
virtually identical. 

Figure 19 shows another test model that tests the bottle and other elements. 
Figures 20 and 21 show time histories of the bottle pressure, the volume pressure, and the 
two flow rates from the “nr” library calculations. Figures 22 and 23 show the same thing 
for the “er” library. The plots are for all purposes identical. 

Figure 24 shows a test model that includes a runtank. Figures 25 and 26 show 
plots of the tank pressure, gas flow rate, and liquid flow rate for computations using both 
libraries. The plots are almost identical. The numbers on the gas flow rate are very 
slightly different, but the change in gas flow rate with time is the same. 

Figure 27 is a test model that considers saturated fluids with mixed-phase states. 
This case is similar to the case considered in Figure 3. Figures 28 and 29 show pressure 
and temperature histories and the two flow rates. These figures can be compared with 
Figures 4 and 5. The results are seen to be similar. The quality in Figure 28 is seen to 
approach 1 and to switch to -1 as the last liquid in the volume boils off. 
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Table 7. Macro Library for nr Library Bottle Module 


C SUBROUTINE BOTLOO ( IPRPL , IUPDAT , MODN , NODI 

C $ NOD2 , VOL , WOUT , HTOUT 

C $ RHOBOT ( UTBOT , HTBOT , DRDT , 

C $ DUDT , TAUCR , TAUCU ) 

C implicit real*8 (a-h, o-z) 

c* ***************************************** **************************** 

C %BEGIN CLASS BOTLOO * 

C ** 

C SUBPROGRAM BOTLOO UNCLASSIFIED * 

C * 

C MECHANICAL ENGINEERING * 

C MISSISSIPPI STATE UNIVERSITY * 

C MISSISSIPPI STATE, MS 39762 * 

C FOR NASA STENNIS SPACE CENTER * 

C %END CLASS BOTLOO * 

c****************** ***************************** *********************** 

C %BEGIN PURPOSE BOTLOO * 

C * 

C ENERGY AND CONTINUITY ANALYSIS OF AN ADIABATIC * 

C BOTTLE WITH ONE EXIT MASS FLOW. * 

C (ADAPTED FROM VOLMOO) * 

C * 

C %END PURPOSE BOTLOO * 

Q * ***************************************** **************************** 

C %BEGIN HISTORY BOTLOO * 

C * 

C WRITTEN 03/17/93 BOB TAYLOR * 

C * 

C * 

C %END HISTORY BOTL00 * 

c************************************************************ ********** 

C %BEGIN SCHEMATIC BOTLOO * 

C * 

c * 


RHOBOT 

UTBOT 

HTBOT 


WOUT NODI 

> HTOUT NOD2 


C %END SCHEMATIC BOTLOO * 

£★*********★*********************************************************** 
C 

c 

stop sort 

call BOTLOO ( nint(PRPB0 — ) , 0 , 'BO—' , 'NODI' , 'NOD2' , 

1 VOLBO — , W1 BO — , HOIBO — , RHOBO — , UT BO — , 

1 HT BO — , DR BO — , DU BO — , TCRBO — , TCUBO — ) 

DR BO— = DR BO — 

DU BO— = DU BO — 

TCRBO — = TCRBO — 

TCUBO— = TCUBO — 

C 

derivative of, RHOBO — = DR BO — 
derivative of, UT BO — = DU BO — 
resume sort 
C 

C Compute Fluid Properties 

C 

C The fluid map indexes are _ 


43 



Table 7. Macro Library for nr Library Bottle Module 


c 

c 

c 

c 

c 

c 

stop 


1 

2 

3 

13 


= H2 
= N2 
= 02 
= He 


sort 

IF ( ICCALC .EQ. 1) THEN 

Use initial values of pressure and temperature to set IC's 
P BO— = PVGBO — 

TTTBO-- = TVGBO — 

VQ BO— = VQGBO — 

IF (VQ BO — . LT . 0 . ODOO ) THEN 
CALL FPROP (NINT (MAPBO — ) , 1, IUNIT, 

1 P BO — , TTTBO — , VQ BO — , d # UT BO — , 

2 HT BO — , S, CP BO—, CV, V, th, 

3 SO, di) 

ELSE IF ( (VQ BO— .GE.O.ODOO) .AND. (VQ BO— . LE . 1 . ODOO) ) THEN 
IP = 1 
IL = -1 

CALL FPROP (NINT (MAPBO — ) , 0, IUNIT, 

1 P BO — , TTTBO—, VQ BO—, d, UT BO—, 

2 HT BO — , S, CP BO — , CV, V, th, 

3 SO, di) 

ELSE 

WRITE ( * , * ) 

WRITE ( * , * ) ' IMPROPER VALUE FOR VQ BO— 7 

STOP 
END IF 

CALL EZSETV (RHOBO — , d/CF00T**3) 

ELSE 

Compute properties based on u and rho 

d = RHOBO— * CF00T**3 

CALL FPROP (NINT (MAPBO — ) , 5, IUNIT, 


1 

2 

3 


P BO—, 
HT BO—, 
SO, di) 


TTTBO — , VQ BO — , d, UT BO — 
S, CP BO — , CV, V, th. 


ENDIF 


P BO — 


P BO — 

VQ BO — 

= 

VQ BO— 

PI1B0 — 

=r 

P BO — 

HT BO — 

= 

HT BO- 

HI1B0 — 

= 

OT BO— 

TTTBO — 

= 

TTTBO — 

TI1B0 — 

= 

TTTBO — 

CP BO — 


CP BO — 

GM BO — 

- 

CP BO— /CV 

GIlBO — 

= 

GM BO — 

SI1B0 — 

= 

S 

RM BO- 

- 

V/cvis 

MI 1B0 — 

= 

RM BO — 

RG BO — 

= 

(CP BO— - CV) 

RG1B0 — 

= 

RG BO — 

RH1B0 — 


RHOBO— 


RJ 


resume sort 
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Table 8. Macro Library for nr Library Inlet Module 


C MACRO source for "IN" - easyrockets exit module 
C Get the map index 
C 1 = para hydrogen , H2 

C 2 = nitrogen , N2 

C 3 = oxygen # 02 

C 13 = Helium , He 

stop sort 
C 

C compute properties as functions of P and T 
C 

VQ IN-- = VQGIN — 

IF (VQ IN— .LT.O.ODOO) THEN 
NISTl = JIDNNT (MAP IN — ) 

NIST2 = 1 

NIST3 = I UNIT 

CALL FPROP (NISTl , NIST2 , NIST3 , PININ— , TTTIN-- , VQ IN—, 

1 d, UININ — , HININ— , S, CP IN—, CV, RM IN—, th, 

2 SO, di) 

ELSE IF ( { VQ IN— .GE.O.ODOO) .AND. (VQ IN-- . LE . 1 . ODOO ) ) THEN 
IP = 1 
IL = -1 

NISTl = JIDNNT (MAPIN— ) 

NIST2 = 0 

NIST3 = I UNIT 

CALL FPROP (NISTl, NIST2 , NIST3 , PININ--, TTTIN— , VQ IN—, 

1 d, UININ—, HININ—, S, CP IN—, CV, RM IN—, th, 

2 SO, di) 

ELSE 

WRITE ( * , * ) 

WRITE (*,*) ' IMPROPER VALUE FOR VQ IN — ' 

STOP 
END IF 


PIlIN — 

= 

PININ— 

VQ IN- 

= 

VQ IN — 

HI 1 IN— 

= 

HININ— 

RHlIN — 


d / CFOOT ** 3 

GI1IN — 

= 

CP IN— /CV 

RG1IN — 

= 

(CP IN CV) *RJ 

MI1IN— 

= 

RM IN— /CVIS 

THIN— 

= 

TTTIN— 

Wl IN— 

= 

W IN— 


resume sort 
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Table 9. Macro Library for nr Library Pipe99 Module 


C SUBROUTINE PIPE99 ( IPRPL , IUPDAT , MODN , NODI 

C $ NOD2 , DIA , ERUF , RKF 

C $ RLEN , PTIN , PTOUT , RHOIN , 

C $ RHOOUT , RMUIN , RMUOUT , W ) 

C implicit real*8 (a-h, o-z) 

c********************************************************************** 

C %BEGIN CLASS PIPE99 * 

C * 

C SUBPROGRAM PIPE99 UNCLASSIFIED * 

C * 

C BOB TAYLOR * 

C STENNIS SPACE CENTER/MISS. STATE UNIV * 

C * 

C %END CLASS PIPE99 * 

Q* ********************************************************************* 

C % BEG IN PURPOSE PIPE99 * 

C * 

C CALCULATE INCOMPRESSIBLE FLUID FLOW THROUGH A PIPE * 

C WITH A LOSS— FLOW COEFFICIENT INTERNALLY COMPUTED. * 

C (ADOPTED FROM PIPE01) 

C * 

C %END PURPOSE PIPE99 * 

c********************************************************************** 
C %BEGIN HISTORY PIPE99 * 

C * 

C WRITTEN 06/09/93 BOB TAYLOR * 

C * 

C %END HISTORY PIPE99 * 

Q********************************************************************** 

C %BEGIN SCHEMATIC PIPE99 * 

C * 

C ================ * 

C NODI > W NOD2 * 

C ================ * 

C * 

C %END SCHEMATIC PIPE99 * 

Q*********************************** **************************** ******* 

stop sort 

RMIPA — = M02PA — 

RMOPA — = Mil PA — 

CALL PIPE99 ( nint ( PRPPA — ) , 0 , 'PA—' , 'NODI' 

$ ' NOD2 ' , DIAPA— , ERFPA — , RKF PA — , 

$ RL PA— , PUPA— , P02PA— , RH1PA— , 

$ RH2PA — , RMIPA— , RMOPA— , W PA— ) 


W1 PA— = W PA— 
W2 PA— = W PA— 
resume sort 
stop sort 

HI 2 PA — = HIlPA — 
HOI PA — = H02PA— 
GI2PA — = GI 1 PA- 
GO 1 PA — = G02PA— 
TI2PA — = Til PA — 
TOIPA — = T02PA — 
resume sort 




Table 9. Macro Library for nr Library Pipe99 Module 



IPRPL 

, IUPDAT 

, MODN 

, NODI 

NOD2 

, NOD3 

, NOD4 

, NOD5 

NOD6 

, VTOT 

, TWALLI 

, HTIN 

HTOUT 

, WIN 

, WOUT 

, pv 

TV 

, HTV 

, RKV 

, TL 

HTL 

, RKL 

, BETAV 

, RNUV 

ALPV 

, BETAL 

, RNUL 

, ALPL 

RHOV 

, WAP 

, UTV 

, UTL 

DRVDT 

RHOL) 

r ****** ** 

, DVVDT 
:********** 

, DUVDT 
**********v 

, DULDT 

t ******* i 


C %BEGIN CLASS RTANK * 
C * 
C SUBPROGRAM RTANK UNCLASSIFIED * 
C * 
C * 
C BOB TAYLOR * 
C STENNIS SPACE CENTER/MISSISSIPPI STATE UNIV. * 
C * 
C %END CLASS RTANK * 

Q* ********************************************************************* 

C %BEGIN PURPOSE RTANK * 
C * 
C THIS MODULE MODELS THE DYNAMICS OF A RUN TANK * 
C PRESSURIZATION WITH VAPOR ULLAGE. INDEPENDANT * 
C MASS AND ENERGY BALANCES FOR THE VAPOR AND LIQUID * 
C ARE USED TO DETERMINE DENSITY AND INTERNAL ENERGY * 
C DERIVATIVES . * 
C * 
C %END PURPOSE RTANK * 


C %BEGIN HISTORY RTANK * 

C * 

C WRITTEN 05/29/96 * 

C * 

C * 

C %END HISTORY RTANK * 

c* ********************************************************************* 

C %BEGIN SCHEMATIC RTANK * 

C * 

C NODI WIN * 

C | < NOD2 HTIN * 

C 1 NOD3 VAPOR * 

C | | PROP3 * 

C VAPOR * 


NOD4 WOUT 
NOD5 HTOUT 
NOD6 LIQ 

PROP5 


LIQUID 


C %END SCHEMATIC RTANK * 

c********************************************************************** 

C 

stop sort 
C 

C VAPOR SECTION 

C 

C The fluid map indexes are 
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Table 9. Macro Library for nr Library Pipe99 Module 


C 2 = N2 

C 3 = 02 

C 13 = He 

C 

IF ( ICCALC . EQ . 1 ) THEN 

CALL FPROP (NINT (MPVRK— ) , 1, IUNIT, PVGRK— , TVGRK — , VQVRK — , 
1 D, UTV, HTVRK—, S, CPVRK— , CV, V, TH; SO, DI) 

PV RK— = PVGRK— 

TV RK— = TVGRK — 

CALL EZSETV(RHVRK— , D/CF00T**3) 

CALL EZSETV (UTVRK — , UTV) 

ELSE 

C CONVERT UNITS ON RHO 

<3 = RHVRK — * CF00T**3 
C 

CALL FPROP (NINT (MPVRK— ) , 5, IUNIT, PV RK— , TV RK— , VQVRK—, 
1 D, UTVRK—, HTVRK—, S, CPVRK—, CV, V, TH, SO, DI) 

ENDIF 
C 

VQVRK— = VQVRK— 

RMVRK — = V/CVIS 

RNVRK — = RMVRK— / RHVRK— 

RKVRK— = TH / CTHK 

ALVRK — = RKVRK— / RHVRK— / CPVRK— 

GMVRK — = CPVRK— / CV 
RGVRK— = (CPVRK— - CV) *RJ 
C 

TBVRK— = TV RK — + 0.001*TV RK — 

CALL FPROP (NINT (MPVRK— ) , 1, IUNIT, PV RK— , TV RK— , VQV1, Dl, 

1 Ul, HI, SI, CPI, CV1, VI, THl , SOI, Dll) 

CALL FPROP (NINT (MPVRK—) , 1, IUNIT, PV RK— , TBVRK—, VQV2 , D2 , 

1 Ul, Hi, SI, CPI, CVl, VI, THl, SOI, Dll) 

BTVRK — = (D2 - Dl)/ (TBVRK— - TV RK— )/(-Dl) 

resume sort 
stop sort 
C 

C LIQUID SECTION 
C 

C The fluid map indexes are 

C 

C 1 = H2 

C 2 = N2 

C 3 = 02 

C 13 = He 

C 

IF ( ICCALC. EQ.l) THEN 

CALL FPROP (NINT (MPLRK— ) , 1, IUNIT, PVGRK — , TLGRK — , VQLRK— , 
1 D, UTL, HTLRK — , S, CPLRK— , CV, V, TH, SO, DI) 

PL RK— = PVGRK— 

TL RK— = TLGRK— 

CALL EZSETV (UTLRK— , UTL) 

CALL EZSETV (WPRK — ,WGRK— ) 

ELSE 

PL RK — = PV RK— 

CALL FPROP (NINT (MPLRK— ) , 8, IUNIT, PL RK— , TL RK— , VQLRK—, 
1 D, UTLRK—, HTLRK—, S, CPLRK—, CV, V, TH, SO, DI) 

ENDIF 
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Table 9. Macro Library for nr Library Pipe99 Module 


RHLRK— = D/CFOOT**3 


VQLRK — 
RMLRK — 
RNLRK — 
RKLRK— 
ALLRK— 
GMLRK— 
RGLRK — 


VQLRK— 

V/CVIS 

RMLRK— / RHLRK— 

TH / CTHK 

RKLRK— / RHLRK— / CPLRK— 
CPLRK — / CV 
( CPLRK- - - CV) *RJ 


TBLRK — = TL RK — - 0.001*TL RK — 

CALL FPROP (NINT (MPLRK— ) , 1 , IUNIT, PL RK— , TL RK— , VQLl, Dl, 

1 Ul, HI , SI, CPI, CV1, VI, TH1, SOI, Dll) 

CALL FPROP (NINT (MPLRK — ), 1, IUNIT, PL RK — , TBLRK — , VQL2 , D2 , 
1 Ul, Hi, SI, CPI, CV1, VI, THl , SOI, Dll) 

BTLRK— = (D2 - Dl)/ (TBLRK— - TL RK— )/(-Dl) 


resume sort 
stop sort 


C 

call RTANK 

1 

1 

1 

1 

1 

1 

1 


( NINT ( PRPRK — ) , 0 
* NOD3 ' , ' NOD4 ' , 
TWIRK— , HI IRK — 
PV RK— , TV RK— 
HTLRK — , RKLRK— 
BTLRK— , RNLRK— 
UTVRK — , UTLRK — 
DULRK — , RHLRK — 


C 


derivative 

of, 

RHVRK — 

= DRVRK— 

derivative 

of, 

WPRK — 

= DWRK— 

derivative 

Of, 

UTVRK — 

= DUVRK-- 

derivative 

of. 

UTLRK- - 

= DULRK— 


C 


resume sort 


stop sort 


, 'NODI' , ' NOD2 ' 

' NOD6 ' , VOLRK— 


, 'RK— ' 

' N0D5 ' , 

, H02RK — 
, HTVRK — 
, BTVRK — 
, ALLRK— 
, DRVRK— 
) 


W1 RK— 
RKVRK— 
RNVRK-- 
RHVRK— 
DWRK-- 


W2 RK— , 
TL RK-- , 
ALVRK— , 
WPRK — , 
DUVRK-- , 


RG1RK— 

= 

RGVRK’ — 

RG2RK — 

= 

RGLRK — 

POIRK — 

- 

PV RK— 

PI2RK — 

= 

PL RK— 

RH1RK — 

- 

RHVRK— 

RH2RK-- 

= 

RHLRK— 

H01RK — 

— 

HTVRK— 

HI2RK — 

= 

HTLRK— 

T01RK — 

— 

TV RK— 

TI2RK — 

= 

TL RK — 

GOIRK — 

~ 

GMVRK— 

GI2RK — 

= 

GMLRK— 

MI1RK-- 

= 

RMVRK — 

MI2RK — 

= 

RMLRK- - 

RK1RK — 

= 

RKVRK — 

RK2RK — 

= 

RKLRK— 


resume sort 




Fig. 14. Input Window for nr-Library Bottle. 
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Fig. 15. Input Window for nr-Library Met Module. 
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Set Uhlts 



Fig. 18. Test Model for Valves and Pipes — nr Library 




Fig. 19. Test Model for Bottle and Volume — nr Library 
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Fig. 27. A Test Model for Saturated Properties — nr Library. 
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SECTION 3 


E-l MODEL DEVELOPMENT UPDATES 

Preliminary EASY/ROCETS dynamic flow simulation models have been 
developed for the H2 and 02 flow systems on the E-l Test Stand (formerly known as the 
Component Test Facility). These models have been updated based on data supplied by 
Mr. Larry deQuay of NASA Stennis Space Center (deQuay, 1997). These new models 
are discussed below. 


High-Pressure Hydrogen 

The high-pressure hydrogen system is shown schematically in Figure 30. High- 
pressure H2 gas is supplied by the bottles. Part of this gas is used in the mixer to temper 
the liquid H2 before it is fed into the test articles. Parallel paths are allowed for flow to 
Test Cells 2 and 3. Figure 31 shows the EASY/ROCETS model for the system prepared 
in the old er-library . Controls are supplied for the runtank pressure and the liquid and 
mixer-gas flow rates. Figure 32 shows the same model developed for the new NIST-code 
nr-library. 

The handwritten notes for the model development are given in Appendix I. 

Figures 33-35 show the results for a trial run using the er-library model. In Figure 
33, the runtank pressure is seen to be held almost constant by the pressure control system. 
Figure 34 shows the flow rates for the runtank pressurant gas, the liquid, and the mixer 
gas. The mixer-gas and liquid flow rates are maintained to set schedules by the 
controllers. Figure 35 shows the control valve flow coefficient for the pressurant gas. 

The two bumps in the curve between 6 and 8 seconds are responses to scheduled bumps 
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in the liquid flow (see Figure 34). Figures 36-38 show the same results for the nr-library. 
As they should be, they are seen to be nearly identical to Figures 33-35. The small 
differences in the starting values are due to the difficulty of setting the two models to 
exactly the same initial conditions. 

High-Pressure Oxygen 

The high-pressure oxygen system is shown schematically in Figure 39. The gas 
bottles feed nitrogen gas to the runtank as a pressurant. The liquid oxygen is supplied to 
Test Cells 1 and 2. Figure 40 shows the EAS Y/ROCETS model for the system using the 
er-library. Figure 41 shows the equivalent model using the nr-library. The handwritten 
notes for the model development are given in Appendix I. 

Figures 42-43 show results of a sample run using the er-library model. Duplicate 
runs using the nr-library model were practically identical. 

Low-Pressure Oxygen 

Figure 45 shows a schematic diagram of the low-pressure oxygen system. Liquid 
02 is supplied in parallel runs to Test Cells 2 and 3. The nitrogen pressurant is regulated 
to 2500 psi upstream of the runtank pressure control valve. Figure 46 shows the 
EAS Y/ROCETS model developed in the er-library. Figure 47 shows the same model 
developed in the nr-library. Figures 48-49 show results of a sample run using the nr- 
library model. Similar runs using the er-library were virtually identical. The handwritten 
notes for the model development are given in Appendix I. 
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Low-Pressure Hydrogen 

The low-pressure hydrogen system is shown schematically in Figure 51. This 
system has not changed since the previous models were developed (Follett and Taylor, 
1996) but models are included here for completeness. Please refer to the cited report for 
a discussion of the model development. Figure 52 shows the model developed using the 
er-library, and Figure 53 shows the same model using the nr-library. Figures 54-56 show 
sample runs using the nr-library model. Equivalent runs using the er-library model were 
practically identical. 


68 



Fig. 30. How Schematic for the High-Pressure Hydrogen System 
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Fig. 32 EASY/ROCETS Model for the HPH2 System Using nr-Library. 
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Fig. 34. Flow Rates versus 
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Fig. 35. Pressure Control Valve Flow Coefficient for the HPH2 System— er-Library. 
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Fig. 38. Pressure Control Valve Flow Coefficient for the HPH2 System — nr-Library. 



Fig. 39. Flow Schematic for the High-Pressure Oxygen System. 
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Fig. 40. EASY/ROCETS Model for HP02 System— er-Library. 
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Fig. 41. EASY/ROCETS Model for HP02 System — nr-Library. 
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Fig. 45. Flow Schematic for the Low-Pressure Oxygen System. 
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Fig. 46. EASY/ROCETS Model for LP02 System— er-Library. 
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Fig. 47. EAS Y/ROCETS Model for LP02 System — nr-Library. 






GENERAL SCHEMATIC FOR THE 
LIQUID HYDROGEN RUN SYSTEM DFFSM 
(LH-DFFSM-1 .1 .8 Apl) 




' NOTES: 

# Any number of bottles can be selected 

<&> Run Tank can be spherical or it can be cylindrical with 
hemispherical/ellipsoid heads 

* Five additional branches can be added for a maximum of 
ten branches; one to ten branches can be selected 

- One to 1 5 line segments in series can be selected for lines 
between nodes shown as solid black dots. Parallel runs can 
be selected for any segment; Labeled nodes have preset \ 
node designations in the DFFSM; other nodes are user * : 
selected, including nodes between user selected line 
segments that are between nodes shown on diagram 

** User can select components to be located between line 
segments as well as sizing/flow parameters for each 
respective component. 

- Gas can be bled off of system at any node (where the user 
selects no component) upstream of the pressure control 
valve 
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Fig. 51. Flow Schematic for the Low-Pressure Hydrogen System. 
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Fig. 52. EAS Y/ROCETS Model for LPH2 System — er-Library. 
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Fig. 55. Pressurant Gas and Liquid H2 Flow Rates LPH2 System — nr-Library. 






SECTION 4 


NEW MODULE DEVELOPMENT 

New models have been developed for an Allen-Bradley SLC 5/03 PID programmable 
logic controller and for a control valve with built in time delay on the actuation. 

Control Systems Enhancements 

The existing EAS Y/ROCETS models for CTF have used general-purpose components to 
implement the PID (proportional plus integral plus derivative) controllers. The actual control is 
planned to implemented using programmable logic controllers (PLC’s). In particular, an Allen- 
Bradley SLC-5/03 PLC has been specified by NASA personnel as a representative device. This 
section details the development and use of an EAS Y/ROCETS module which provides a model 
of the PID control function of the AB SLC 5/03, which actually is constant over the SLC 5 series 
of PLC’s. 

Derivation 

A PID controller generally implements some discrete approximation to the equation 

dE /i \ 

Edt + K d —— 
o D dt 

Implementation on a digital computer such as the PLC requires a numerical approximation to this 
exact equation. In addition, the derivative term must be filtered to limit the effects of high 
frequency noise terms. A commonly used scheme to approach this problem is to incorporate a 
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low-pass filter on the derivative term [Phillips and Harbor, 1996]. If Equation (1) is converted to 
the Laplace domain, we get 

CV(s)= K p E(s)+K I jE(s)+K D sE(s) (2) 


where the term (1/s) represents an integration and the term (s) represents a differentiation effect. 
Applying the low-pass filter to Equation (2) yields 


CV(s) = K p E(s) +K,j E(s) + K d 


1 + 


■E(s) 


co f 


(3) 


The variable co f represents the cutoff frequency of the low-pass filter, which has transfer function 


G/0)= \/ (4) 

1+ V 

/ COf 

The choice of co f is dependent on the desired frequency response of the PID. If co f is chosen too 
small, the derivative term will be ineffective due to the rolloff of the filter G f (s). If it is chosen 
too large, high frequency noise will cause problems due to the amplification effect of the 
derivative term. 

The Allen-Bradley SLC 5 series of PLC’s implements a filtered form of the PID equation 
as described in Equation (3). However, the details fo the implementation, including the 
placement of the derivative term filter cutoff frequency, are considered proprietary information 
by Allen-Bradley personnel. To date, no clear information regarding the details of the 
implementation have been forthcoming. 

In order to provide a component for use in the EAS Y/ROCETS environment, certain 
assumptions have been made regarding the detailed operation of the PLC version of the PID 
controller. First, it is assumed that the calculations are simply discretized state models of the 
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integral and derivative terms of the control variable as shown in Equation (3). Additionally, the 


filter frequency co f is placed at the frequency of operation of the PID, which is specified as the 
scan time of the ladder diagram. Note that this is a reasonable assumption based on sampling 
theory. It is well-known that the sample rate should be twice the highest frequency of interest. 

By filtering above this, we tend to limit the high-frequency terms that are considered to be noise. 
It should also be noted that this choice was made after observation of the operation of an example 
system with a variety of scan times and gain constants. 

By splitting Equation (3) into three pieces, we can calculate the proportional term, the 
integral term, and the derivative term individually, then combine them to give the complete 
answer for the control variable. The three terms are given in Equation (5): 

CV p (s ) = K p E(s ) 

CV,(s) = — E(s) (5) 

£ 


CV D (s) = 


K D co f 

s+cOf 


E(s ) 


Writing continuous state equations from these terms, using the integral and derivative terms CY! 


and CV D as states, we obtain: 


CVj = KjE 

J ' ( 6 ) 

CV D = -0) f CV D - K D a> f 2 E 

with 

CV = CVj + CV D + K P E (7) 

Equations (6) show the state equations and Equation (7) gives the output equation. Next, the 
continuous equations are discretized, yielding 
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( 8 ) 


CVj (k + 1) = [1] CV t (k) + [K, T]E(k) 

CV D (^ + l) = [^ r ]Cy D (A:) + (-^ / ^)(l-^ T )£(A:) 

Note that the sampling period for the discretization is represented as t, and corresponds to the 
scan time of the PID controller. Using the previously mentioned assumption that co f =l/x, we get 

CV J (k + 1) = [1 ]CVj ( k ) + [K I t]E(k) 

CV D (k + 1) = [e~'] CV D (k) + (-K„/ m-e- l )E(k) <9) 

The output equation for the system has the same form as the continuous output equation, so that 

CV(Jfc) = CVj (k) + CV D (k ) + K p E(k) (10) 

Equations (9) and (10) give the discrete-time state-space description of the operation of a 

standard PID controller. Customization of the PID for the AB-SLC 5 series of PLC’s is 

accomplished by means of a conversion to a set of dependent gains, which provide an alternate 

description of a PID controller. For the SLC 5, the three terms which specify the operation of the 

PID controller are K c (proportional gain term), Tj (integral gain or reset term), and T d (derivative 

gain or rate term). Conversion between the two sets of gains is as follows: 


K P = K C 

K= 

1 607: 

K d = 60 K c T d 


( 11 ) 


where Tj has units of repeats/minute and T d has units of minutes. Additionally, the values are 
entered into the SLC 5 as integer values only, so a range and gain enhancement specification 
allows K^. and T ; to be entered using a multiplier of either 10 or 100. T d , however, is always 
entered using a multiplier of 100. The conversions can be represented using the variable “gain r 
to show the adjustable gain (either 10 or 100). 
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(12) 


K P = 


K,= 


K d = 


K_ 

gain 

A_ 

607: 

( k„ Y z 


y gain J 


100 


1(60) = 0.6^- 
I gain 


So, by using the conversions in Equation (12), we can implement the state and output equations 
of Equations (9) and (10). Component AB of the EASY/ROCETS library er is an 
implementation of these equations. Figure 57 shows the icon for the AB component, and Figure 
58 shows a portion of the Input Specification Table. 


AB-SLC-5/03 PID 


AB 

C = 1 10 
Ti= 367 
Td=2 


Fig. 57. Icon for Module AB (Allen-Bradley SLC-5/03) 
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Fig. 58. Component Data Table for Module AB 






EASY/ROCETS Usage of AB 


The principal design objective for the AB module was to develop a model that has I/O 
characteristics that are as close as is possible to those of the actual device. With that in mind, the 
inputs for the AB module are as shown in Table 10. 


Variable 

Description 

Units 

Default Value 

SP 

setpoint 

Engineering units or scaled for PID 

1.0 

RG 

reset and gain range 

unitless 

1.0 


enhancement bit 



C 

controller gain 

unitless (/10 or/100) 

.99999 

TI 

reset term 

minutes/reset (/10 or /100) 

.99999 

TD 

rate term 

minutes (/100) 

.99999 

TAU 

loop update time 

seconds 

1.0 

PV 

process variable 

Engineering units or scaled for PID 

.99999 

SC 

scale setpoint flag 

unitless 

0 

AP 

scaling rate 

unitless 

1.0 

APO 

scaling offset 

unitless 

1.0 

BI 

bias 

bias or feedforward term 

0 

DB 

deadband 

Engineering units or scaled for PID 

0 


Table 10. Inputs for AB PID Module in EASY/ROCETS 
The variables defined in Table 1 are given default values as shown in the table. Also, the 
units for each variable are shown where appropriate. For SP, PV, and DB, the units can either be 
in engineering units, or in units scaled for the PID calculations. It should be noted that whether 
the variables are in scaled form or engineering units, the same numerical accuracy applies. 

Where default values do not make sense, the value .99999 is used. Note that if SC is set to 0, the 
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values of the scaling constants AP and APO do not matter. The primary output of the module is 
CV, the control variable output state. Other output states provided include the current error (ER), 
the control variable contributions due to the integral and derivative terms (CVI and CVD, 
respectively), and the PHD versions of the process variable and the control variable (PVP and 
CVP, respectively). In addition, there is a state ACV that is the control variable without bias. It 
is an intermediate variable used in the calculation of the eventual output CV. 

Demonstration of AB Usage 

To demonstrate the use of the new AB module, a simple system is shown in Figure 59. 
The system is a third-order transfer function with poles at 0, -1, and -2. It is desired to maintain a 
constant setpoint of 1, using PID control. This system has been analyzed using the actual PLC 
hardware as well as in the EAS Y/ROCETS environment. For the hardware setup, the dynamic 
system was simulated in real-time using a personal computer running VisSim/32 and an I/O card 
to provide real-time input and output. 

The SLC-5/03 which was provided by NASA personnel for testing uses a current loop 
output module, which requires a conversion to a voltage signal in order to work properly with the 
VisSim simulation of the system plant. A 500 Q resistance was used, which produced a 10V 
output for the maximum output of 20 mA. In order to obtain a bipolar output, a shift of -5V was 
included inside the VisSim simulation. This setup is shown in Figure 60. Note that the control 
variable and the process variable are shown in their various configurations as CV and PV. The 
simulation diagram of the EAS Y/ROCETS version of the system is shown in Figure 61. 
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Fig. 59. Basic Block Diagram of a Demonstration PID Control Loop 



(0-20mA) 


Fig. 60 PID Control Loop Simulation with SLC-5/03 Hardware in the Loop 





Fig. 61. EASY/ROCETS Block Diagram for Demonstration System 




In order to compare the outputs of the system with the SLC-5 hardware in the loop and 
the system which is entirely simulated in the EAS Y/ROCETS environment, the following values 
were used in both systems: 


SP = 819.2 

(PID form of setpoint corresponding to 1 V) 

RG = 1 

(use the extended range, which divides C and TI by 100) 

C= 110 

(K c = 1.1, then multiply by 100) 

TI = 367 

(Tj = 3.67, then multiply by 100) 

TD = 2 

(T d = 0.02, then multiply by 100) 

TAU = 0.01 

(scan time of 10 msec.) 

SC= 1 

(scale the setpoint and the deadband) 

AP = 819.2 

(scaling rate) 

APO = 0 

(scaling offset) 

BI = 8192 

(bias the output by Vi full scale) 

DB = 0 

(no deadband for this example) 

Figure 62 shows the output of the hardware in the loop version of the system with 


results of the EAS Y/ROCETS simulation. There are obviously some differences in the two 
responses. These differences are shown in Table 1 1, along with the relative error for each 
parameter. These differences are attributable to the lack of a detailed knowledge of the actual 
equation used in calculation of the PID control variable. The detailed algorithm used in the 
hardware would allow precise calculation of the responses, including quantization effects 
associated with the integer arithmetic employed. However, as has been mentioned previously, 
the details of these calculations are proprietary information, and are not available at this time. 
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Output of the Plant 



Fig. 62. Demonstration System Response for Hardware Testing and 
EASY/ROCETS Simulation (No Deadband) 
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Parameter 

Hardware Value 

Simulated Value 

% Error 

Time-to-Peak 

3.220 

2.578 

19.9% 

Peak Value 

1.254 

1.195 

4.7% 

Steady-State Value 

1.022 

1.003 

1.9% 


Table 11. Comparison of the Hardware and Simulated Responses for No Deadband. 

Next, the system was simulated with exactly the same values as before, except that a 
deadband was inserted with magnitude 100. The results are shown in Figure 63, and the results 
are summarized in Table 12. The responses again are similar, but have some significant 
differences. It is believed that the hardware response has more of an oscillatory response due to 
the fact that the peak value was higher and the deadband played a more significant role in the 
response. If the amplitude scaling is adjusted so that there is a smaller peak value, the response 
does not go as far out of the deadband, and does not develop the larger control variable that 
results from the response shown. 


Parameter 

Hardware Value 

Simulated Value 

% Error 

Time-to-Peak 

3.125 

2.578 

17.5% 

Peak Value 

1.296 

1.157 

10.7% 

Steady-State Value 

1.033 

1.004 

2.8% 


Table 12. Comparison of the Hardware and Simulated Responses with Deadband 
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g. 63 


Demonstration System Response for Hardware Testing and 
EASY/ROCETS Simulation (Deadband of 100) 





Conclusions 


It is evident from the previous simulations that the EASY/ROCETS module AB gives a 
reasonable approximation to the response of the actual Allen-Bradley SLC-5/03 hardware. The 
differences are accounted for by the lack of detailed knowledge of the hardware implementation 
of the PID control law. However, the general form of the response is reproduced by the 
simulation, and can still provide simulations which can be used in the design of PID control laws 
for use in systems which can be modeled by EASY/ROCETS. Once the PID controller is 
designed and simulated in EASY/ROCETS, the results can be used as a baseline from which to 
fine-tune the details of the controller. This will significantly reduce the number of iterations 
required to choose the controller constants, with resulting savings in both time and cost. 
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Valve Module with Delay f Component XD) 

The valve module with delay is adapted from one developed by Steve Poulton to model 
compressible fluid flow through a valve. The valve actuation has been modified to include a 
delay term to give a more realistic model of response times. The following section shows a brief 
derivation of the code and discusses the usage of the module in the EASY/ROCETS environment 


Derivation and Use 

A preliminary analysis of the time delays associated with changes in valve stem position 
has been made by Follett and Taylor (1996). The delay time was modeled by a lag function, 
which gives an exponentially varying response time in response to a step-type input. Figure 64 
shows a schematic of this implementation, and Figure 65 shows response curves comparing the 
command given to the valve and the simulated response of the valve. Note that in the model 
shown, the fully open flow coefficient is the parameter being varied rather than the percent open 
coefficient. By connecting the control input to the variable representing the percent open, the 
valvestem position could be directly controlled. As is seen in the resulting plot, a much more 
realistic response curve is obtained by this technique. The variable which controls the speed of 
the response is the time coefficient in the denominator of the lag function. 

The details of the operation of this new module are clearly evident when the multiplication 
by the lag function is considered. A lag function has a general form of 



1 

ts+1 


(13) 


and can be represented in the time domain by a multiplication by a time delay of 


l delay 


= e 


-tlx 


(14) 
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Fig. 64. Block Diagram Setup for Demonstration Case for Valve Delay Module 
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with the variable x representing the time constant of the lag. In a period of time equal to four time 
constants, the response has traveled through 98% of the change that it will make. So, a valve with 
this lag function built in will change by 98% of the commanded change within this time period. 

The new valve module accepts as an input the maximum actuation time specification for 
the valve. Internally, this number is divided by four to obtain the time constant to be used in the 
module. Additionally, the default value of zero can be used to obtain ideal, or instantaneous 
responses. Also, the user can specify whether the delay is to be applied to the fully open flow 
coefficient or the valve percent open. Either one of these parameters can actually be used as the 
control variable for the valve, and therefore must be able to be delayed by the valve actuation 
time. A copy of the valve icon for the new delay valve is given in Figure 66. 


Valve 98 with time delay 



Fig. 66. EASY/ROCETS Component Icon for Valve with Delay 
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SECTION 5 


CONCLUSIONS AND RECOMMENDED FUTURE WORK 

Several advances in the state of the EASY/ROCETS dynamic fluid flow 
simulation package have been presented in this report. Incorporating the complete NIST- 
12 property database into the EASY/ROCETS library has significantly advanced the 
thermodynamic property capabilities of the system. This involves a complete new library 
“nr-library,” which has been redesigned for more natural data entry for the initial 
conditions. New modules have been developed for the Allen-Bradley SLC 5/03 
programmable logic controller PID system and for a valve with time delay actuation. 
Models of the E-l Test Stand subsystems have been updated to the most recent design 
data. Using the new nr-library and the revised er-library, these updated models now work 
for all pressure levels. 

EASY/ROCETS has reached a level of maturity such that it is useful for 
performing simulations of flow systems in support of facility design, test scenario 
development, facility modification and facility control. Further augmentation of the 
package should incorporate feedback from system engineers at SSC. The future 
usefulness of EASY/ROCETS depends to a great extent on the perceptions of the 
individual engineers regarding ease of use and validation of the output. Once system 
models are used on a regular basis, clear plans can be made for specific upgrades or 
changes to the software. 

Some specific recommendations for future work can be made at this time. 
EASY/ROCETS should be worked heavily to develop control system designs and to 
work on the initial tuning of the PID controllers. The package has not been used in this 
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mode much to date and both the present researchers (Follett and Taylor) and our 
colleagues at SSC need to exercise the package in this mode to determine how modeling 
can help with the control system and what modeling strategies need to be worked out, to 
tune controllers for example. Also, more systems at SSC should be modeled in this 
system. 

Validations should still be performed by comparison with actual facility data 
accumulated during test runs. As tests are made on systems that have been modeled, the 
actual results should regularly be compared with the model predictions, and the 
comparison used to increase the fidelity of the models. Detailed analysis by simulation 
experts should occur on a regular basis. 
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APPENDIX I 


HANDWRITTEN NOTES ON MODEL DEVELOPMENT 

HPH2 119 

HP02 135 

LP02 I 43 
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High-Pressure Hydrogen System 
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